#include "FiniteDiff/MixedDer.h"

#if DIM == 2

template <class T>
void MixedDer::apply(Real alpha, const Tensor<Real, Dim> &phi, const T &rhs, Tensor<Real, Dim> &LofPhi) const
{
  auto dx = rd.spacing();
  loop_box_2(rd, i, j) {
    Real xx = (-phi(i-2, j) + 16*phi(i-1, j) - 30*phi(i, j) + 16*phi(i+1, j) - phi(i+2, j)) / (12*dx[0]*dx[0]);
    Real yy = (-phi(i, j-2) + 16*phi(i, j-1) - 30*phi(i, j) + 16*phi(i, j+1) - phi(i, j+2)) / (12*dx[1]*dx[1]);
    Real xy = 0.0;
    Real c[] {1.0, -8.0, 8.0, -1.0};
    int  d[] {-2, -1, 1, 2};
    for(int l = 0; l < 4; ++l)
      for(int k = 0; k < 4; ++k)
        xy += phi(i+d[k], j+d[l]) * c[k] * c[l];
    xy /= (12.0 * 12.0 * dx[0] * dx[1]);
    LofPhi(i, j) = alpha * (A[0]*xx + A[1]*xy + A[2]*yy) + rhs(i, j);
  }
}

void MixedDer::apply(const Tensor<Real, Dim> &phi, Tensor<Real, Dim> &LofPhi) const
{
  apply(1.0, phi, [](int, int) { return 0.0; }, LofPhi);
}

void MixedDer::apply(const Tensor<Real, Dim> &phi, const Tensor<Real, Dim> &rhs, Tensor<Real, Dim> &rsd) const
{
  apply(-1.0, phi, rhs, rsd);
}

void MixedDer::relax(const Tensor<Real, Dim> &phi, const Tensor<Real, Dim> &rhs, Tensor<Real, Dim> &JofPhi, Real w) const
{
  auto dx = rd.spacing();
  loop_box_2(rd, i, j) {
      Real xx = (-phi(i-2, j) + 16*phi(i-1, j) + 16*phi(i+1, j) - phi(i+2, j)) / (12*dx[0]*dx[0]);
      Real yy = (-phi(i, j-2) + 16*phi(i, j-1) + 16*phi(i, j+1) - phi(i, j+2)) / (12*dx[1]*dx[1]);
      Real xy = 0.0;
      Real c[] {1.0, -8.0, 8.0, -1.0};
      int  d[] {-2, -1, 1, 2};
      for(int l = 0; l < 4; ++l)
        for(int k = 0; k < 4; ++k)
          xy += phi(i+d[k], j+d[l]) * c[k] * c[l];
      xy /= (12.0 * 12.0 * dx[0] * dx[1]);
      JofPhi(i, j) = (A[0]*xx + A[1]*xy + A[2]*yy - rhs(i, j)) / (30.0/12 * (A[0]/(dx[0]*dx[0]) + A[2]/(dx[1]*dx[1])));
      JofPhi(i, j) = JofPhi(i, j) * w + phi(i, j) * (1-w);
    }
}

#endif // DIM == 2
